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Abstract 

In this paper a non-iterative method for the numerical generation 
of orthogonal curvilinear coordinates for plane annular regions between 
two arbitrary smooth closed curves has been developed. The basic 
generating equation is the Gaussian equation for an Euclidean space 
which has been solved analytically. The method has been applied in 
many cases and these test results demonstrate that the proposed method 
can be readily applied to a wide variety of problems. The method can 
also be used for simply connected regions only by obtaining the solution 
of the linear equation (19) under the changed boundary conditions. 

Details on the work reported in this paper are available in Reference [1] . 
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Fundamental Ideas of the Method 

All methods of numerical coordinate generation in a two-dimensional 
plane and classified under the method of "elliptic equations" (Refs. 

[2]- [10]), have depended invariably on the solution of Poisson equations 
for the curvilinear coordinates £(x,y) and n(x, y) : 


g (g li r 22 2g 12 r i2 + g 22 r il ) 


6 22 


v 2 n = 


-=■ (e r 2 

g vg ll 22 


2g 12 r i2 + g 22 r il ) 
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inhere P(£,n), Q(£,n) are arbitrarily specified control functions, the g_j 
are the fundamental metric coefficients, the V\ are the Christoff el 
symbols of the second kind 
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Implicitly equation (1) implies two things: (i) that the coordinates 

for the same domain can also be obtained by solving the Laplace equations 

V 2 £ = 0, v 2 n = o (3) 
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and (ii) since the have first partial derivatives of g_ in them, 

equation (1) can also be interpreted as providing a set of constraints 
or relations among the g and their first partial derivatives. 

In this paper we present another method based on elliptic 
equations and state the problem as follows. 

The three functions g^, £^2 ’ §22 curv ^^-^ near coordinates 

i 9 r\ define an element of length ds in a plane if the Gaussian 
equation with zero curvature 

/g r 2 /g r 2 

— ( e n ) - a ( . g h2 = o ( , 

an ^ g n ; H K g n } 

holds for every point in the plane, and then the Cartesian coordinates 


are given as 


x = x(£,n) , y = y(C,n) 


Equation (4) is identically satisfied by a function a(£,n) defined as 


-4 


a = — & Ff , a = — “■ rf 0 
4 8 11 11 n 8 11 12 


Specifically, a is the angle of inclination with respect to the x-axis 
of the tangent to the coordinate line n = const, directed in the sense 
of increasing values of the parameter £ . The first partial derivatives 
of x and y are 


x^ - /g 1]L cos a, y - -/g^ sin a 


— — (g 10 cos a + Jg sin a) , y « - — - — (/g cos a - g. 0 sin a) I (5) 
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Then 


cos a dC + — (g^ 2 COS a + sin a )dn] 
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y = -/ [ /g^ sin a - ' (/g cos a - ? sin a)dn] 
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a - -/ ij (r L d? + r i 2 dri) 


( 6 ) 


The inverse relations of (5) are 


£ x * (>^g COS a - g 12 sin a)//gg^ 


C y = -(g 12 cos ct + /g sin a)//gg^ 


n x = sin a 


n y = ’ / S 11 /g cos a 


(7) 


For the case of orthogonal coordinates, the coefficient 


g l2 - 0, i.e.. 


8 12 ' Vn + y c y n ’ 0 


( 8 ) 


which is satisfied by the equations 


x 

n 



y 


n 



(9) 


where F > 0 is a continuous function of g and [11] • 

Referring to Figure 1, let the boundary T ^ of a bounded region 
in an Euclidean two-dimensional space be a simple curve x « x^CO, 
y * y^CC), with a uniformaly turning tangent* In the region ft, let 
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be an annular subregion bounded by the inner boundary and the 
outer boundary The region is to be mapped onto a rectangular 

region R in the £q-plane by a transformation so as to have 




x = x(£,n) 

y = y(£,n) 


n e ± n ± ^ 


where q rt and n are the actual parametric values associated with the 

p CO 

boundaries and respectively, and x,y are periodic in the ^-argument. 
Substituting g^ = 0 in the fundamental equation (4) we get 
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( 10 ) 


where 

8 22 = FZg ll | 

g = (Fg n ) 2 ) (ID 

Before we solve the problem of orthogonal coordinate generation 
based on the elliptic equation (10), we digress and state the following 
results: Following Potter and Tuttle [ 6 ] we assume that the £ -curves 
in the physical xy-plane are free from sources and sinks. This 
condition establishes a unique correspondence between the 5-points on 
each pair of q - const, lines. In the absence of sources and sinks, 
we have 

div[grad ip(q)] = 0 (12) 

where ip(q) is an arbitrary differentiable function of q and grad iKn) 
is oriented along the normal to the curve q = const. Carrying out the 
differential operation in (12) and using the expressions 


| grad n| = 1/^22 




V 


524 


in (12) , we obtain 



3n 


Un > / g^7g^’) = 


d 2 ip , d^ 

dn 2 dn 


Writing = l/v(n) and denoting the arbitrary function of £ due to 
integration by Jlny(C), we obtain the result 


' /g n/ g 22 = V ^)v(n) 

= 1/F (13) 

This result shows that for the case of orthogonal coordinates the 
ratio g^/%22 is a P rociuct of the positive functions p (O and v(n). 

The result in (13) also provides the condition for the two distinct 
families of orthogonal curves 

£ = const., n ° const. 

to divide the physical plane in inf initensimal squares. (See Cohen [12]). 
We now introduce new coordinates (£) and n T (n) as 

S' - JxtOdt, n’ ■ I ^ U4) 

Thus 

S 11 = 8 ll y,y2 ’ g 22 = g 22 v2 


so that 




Defining 


P' = In g 
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it can be shown that 
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(15) 


Using (15) in (10) , we get a much simpler equation 

^ + ^ p - - 0 (16) 

95’ 2 9 n ' 2 

Another important result can be obtained based on (13) . Using 
the orthogonality condition g^ = 0 in (7) , we have 

n = 5 /f, n * - £ /F 
y x x y 


so that 


-k <V F) + i V F) ‘ 0 


(li) 


Carrying out the transformation (14) in (17), we get 


V 2 £ T = 0 


(18) 


Equation (18) provides the uniqueness condition for the solution of 
equation (16) . 

Based on the preceeding analysis we can state that if an exact 
analytic solution of equation (10) can be obtained for F = 1, i.e., 

^22 " then the solution for any other coordinate system £ and r \ , 

where £ = <£ (£) and n = f(n) 5 can simply be obtained by the substitution 
of <j> and f in place £ and r\ respectively. With this scheme in mind, 
we solve the equation 


9 2 P 

H 2 


+ 0 
3n 2 


( 19 ) 


where 


H' " f!' ITT? RP’SKIT^f S f. '\ r ; 


( 20 ) 
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P = Jin g 


under the boundary conditions 


P = P^(0 at n = 0+ 


= P (£) at n = n 


11 


0 < £ c 2 tt 


( 21 ) 


where the subscripts g and <» denote the inner and outer boundaries, 
respectively. The periodicity requirement is that 


P(5,n) = p(^+2tt ,n) (22) 

Further, the ^-coordinate must be such that the equation 

V 2 s = 0 (23) 

is always satisfied. 

A general analytic solution of equation (19) under the boundary 
conditions (21) and the periodicity condition (22) is 


where 


and 


P(5,n) = a + nK + 1 - sinh n(n -n) 
o n=l 00 

(a cos n£ + b sin nF)/sinh nr] 
n n 00 

00 

+ l 1 sinh nr|(c cos n£ + d sin nO/sinh nn (24) 
n=l n v n 00 


K = (c -a )/n 
o o 


oc 


(25) 


2tt 


, _L J 2 ’ % - £ J 0 P e (5)d« 


a o 2rr ' o 


i*There is no loss of generality in setting the parametric value ng = 0. 
The value n must be interpreted as the difference between the actual 
values of r\ at the outer and inner boundaries. 
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1 t ' P.(£) COS n£ d? , b = 
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For orthogonal Coordinates 
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therefore for 
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and consequently 

cosh nCn^-n) 

a(C,n) * oi(£,o) + E , 9 ~ 7 ~T ~T ~ (b cos n£ 

n=l l sinn n n n 


a sin nO 
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, oo cosh nq , v 

+ L t: r-r (c sin n£ - d cos nu 

n=l 2 sxnh n n n n 


cosh n n 

oo oo 

- Z . - — r — (b cos n£ - a sin n£) 

n=l 2 sxnh nv n n 


- Z - -r : — r (c sin n£ - d cos n£) 

n=l 2 sinh n q n n 


Having determined g ^ and a, we can find the Cartesian coordinates 


:(£,n) = x(C,o) + /g 22 sin ct dn 


y(S,n) - y(5,o) + cos a dn 


The preceding solution is for the case when g 22 = i.e., 

F - 1. However, as stated earlier, the solution for any other coordinate 


system £,q in which g 22 ^ can be obtained by replacing £ and n 

in (24), (27) and (28) by <j>(£) and f(q), respectively. This feature can 
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be used to redistribute the coordinate lines in the desired regions. 

Since the functions <f> and f are at our disposal, they play the same role 

as P and Q in equation (1). Further, since the Fourier coefficients 

(26) are invariant to a coordinate transformation £ = <)>(£)> where 

d) (£ ) = 0, <j> (£ ) = 2 tt and £ corresponds to £ = 0, 5 corresponds to 
o m o m 

£ = 2tt, these coefficients need not be recalculated. 

The procedure of transformation from £,n to £,n is as follows. 

On transformation from (£,p) to (£, n), the covariant metric 
coefficients transform through the equation 


- _ jbc^ jhd 

^ " §k£ 3 k 1 ^ 

so that, on using the relations g ^ = g^^ and g = 0, we have 


’ll 


= [(% 2 + (% 2 ]g 


35 


35 


11 


■ 35\ 2 , , 1 Hn 2 


§22 = [(■ _) + <T=) ]gii 

Sn 3n 


(29) 


We now introduce the transformation 

5 = <t>(5) | 

n = f(n) ) (30) 

where the functions $ and f are continuously differentiable and satisfy 
the conditions 

<K5 0 ) = 0, f (rig) = 0 


where £ = 0 corresponds to £ = and q = 0 corresponds to q = n^. 


Defining 


A = ^ , 6 = ^ 
d5 dn 
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we obtain from (30) 


0 2 

g 22 (S,n) = — g l;L (C,n) 


Comparing with (13), we find p * A, V = — . 


(31) 


The salient feature of the preceding analysis is that the solution 
for the case g ^ can be used to obtain the solution for the case 

^22 ^ Sfi by coordinate transformation. 


Before solving any specific problem, it is important first to 
establish an orthogonal correspondence between unique points of the 
inner and outer boundary curves. This condition is satisfied if we 
choose the ^-curves satisfying the equation 

V 2 5 = 0 (32) 

The inner and outer boundary curves are available either in tabular 
or functional form as 

y e = y( V’ y “ = y(x » } (33) 

For equation (32) to be satisfied, we can take £ as the angle traced 
out in a clockwise sense by the common radius of the concentric circles 
in a conformal representation of the inner and outer boundary curves. 

If a and A, respectively, are the radii of the inner and outer circles 
in the transformed conformal plane, then 


n 

i 2 " 

J o 

lx B <{) 

cos £ 

- y 6 «) 
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€]d£ ' 
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2tt 
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cos £ 

- y JO 

sin 

£]d£ 


(34) 
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As is well known, the preceding scheme is an iterative numerical scheme. 
In lieu of this, we have developed a method which is fast and direct, 
and is equivalent to satisfying equation (32), 

We circumscribe circles around the inner and outer boundary curves. 
Two cases arise depending on whether the circles are concentric or 
nonconcentric . 

Case I : If the circumscribed circles are concentric (Fig. 2(a)), then we 

select those values of the ordinates which correspond to the abscissae 

x _ = r cos x = r cos £ (35) 

p s 00 L 

where r and r are the radii of the circumscribed circles in the 
s L 

physical plane. 

Case II : If the circumscribed circles are nonconcentric (Fig. 2(b)), 

then we first use the formula for the conformal transformation of non- 
concentric to concentric circles, Kober [13] , and choose the ordinates 
corresponding to abscissae given by the formula 

x(£) = [ (1 - cy cos £) {x (1 - cy cos £) + cyy sin £ 

t j-i 

+ r (c cos if) - y cos(£-i|0)} 

L 

- cy sin £{y (1 - cy cos £) - cyx sin £ 

Li Li 

- r (c sin ip + y sin(£-ijO)}] 

/ (1 - 2cy cos £ + c 2 y 2 ) (36) 

r ,r * radii of outer and inner circumscribed circles 
L s 

(x_ , y_ ) and (x , y ) = coordinates of the centers 
Li Li s s 


where 



(a) Concentric circumscribed circles Cj 
and C£ of radii r s and r^, respectively, 
with center at the origin. 



(b) Nonconcentric circumscribed circles Cj 
and C 2 of radii r s and rj and centers 
at z s and respectively. 


Figure 2.- Circumscribed circles. 





d2 = (x s' X L )2 + W 2 


y s _y L 

'!>=*- tan-1 C ^~ ) 
s L 


c = [ (d 2 + r 2 - r 2 ) + { (d 2 + r 2 - r 2 ) 2 - 4d 2 r 2 /2r T d 

L S 


L s' 


y = 1 for the outer boundary 


Y = ^° r t ^ le ^ nner ^ oun dary 

s 


Having determined the appropriate sets (x (£) , y (5)) and 

p p 

(x^Q, y (?))> we use (34) to obtain the values of a and A. The 
parametric difference n is connected in some manner with the "modulus’ 1 
of the domain which, however, by itself is a separate problem (see 
Burbea [14 ] and Gaier [15 ] ) . In this work we have defined n based 

oo 

on the knowledge of a and A as discussed above by the formula 

n = &n(-) (37) 

oo a. 


For Figures 3 to 8, we have used the following functional forms 
of (f> and f : 


4>U) = 


2irU-£ 0 ) 


n.(n-y k^V 

f(n) = " n — — — — — zr~ 

n -n e „(n -ti ) 
«>3 K 'o 3 
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so that 


A = 


2tt 

t-? n 

m o 


e = 


n (n-nJ 

^eo K ^ 

— - [l + (n-n g )toK]- -- r - 


n -n 


3 


K P 


where K > 1 is an arbitrary constant, and £ = £ , n = correspond, 

respectively, to £ = 2tt and r) = r^. We treat £ and n as integers so 

that £ = 1, ^ = I MAX, n = 1 , and n = JMAX. Since n is known 
o m 3 00 00 

from (37), hence by specifying the numerical values to K and JMAX 
we can create the desired mesh control in the direction of n* The 
value of K between 1.05 and 1.1 is quite sufficient [16] to have a 
fine grid near the inner boundary. 

The number of terms to be retained in the series (24) is usually 
small for convex inner and outer boundaries, though we have retained 
(IMAX-l)/2 number of coefficients in each computation. This number 
is the optimum number of terms in a discrete Fourier series [ 17 ] 
having IMAX number of points in one period. 

Figure 3 shows the classic case of confocal ellipses with coordinate 
contraction in q. The value of K is 1.05. The orthogonal correspondence 
between 5 -points of the inner and outer boundary has been established 
by using Case I, Eq. (35). 


u 



Figure 3.- Confocal ellipses. Semimajor axes 1.48, 5.0, 
and semiminor axes 0.5, 4.802, respectively. Only 
38 n = Const, lines shown for detail. 
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Figure 4 . - 


A blunt body section with elliptical outer boundary. 
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Orthogonal coordinates for a Joukowsky’s airfoil with slightly rounded 


trailing edge are shown in Figure 6. Eq . (35) is used for orthogonal 
correspondence. Here K = 1.02. 



Figure 6.- Joukowsky's airfoil with slightly 
rounded trailing edge . 


538 


Figure 7 presents orthogonal coordinates for nonconcentric ellipses. 
Centers of the inner and outer ellipses are at (0,0) and (1,0), respectively. 
Here K = 1.01. For orthogonal correspondence Eq. (36) has been used. 
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Figure 8.- Generated coordinates for body having convex, 
concave and straight portions. Placement of outer 
boundary is decided by the radius of the osculating 
circles of the concave portions. 
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Summary of Numerical Experimentation 


In the course of this investigation a number of cases of inner 
and outer boundary shapes and orientations have been tested through 
the developed computer program. The main conclusions are listed below: 

(i) The method works very effectively for smooth and convex boundaries 
of any shape and orientation. 

(ii) For concave boundaries a method similar to that of Eiseman has 

to be used in the placement of the outer boundary to avoid intersecting 
normals. Another remedy would be to introduce some type of 
attraction near the outer boundary facing the concave side of the 
inner boundary. 

(iii) Sharp turns and corners are not admissible and have to be rounded 
to avoid singularities in the metric data. 
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